<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Polynomial discrimination</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/cvxbook/Ch08_geometric_probs/html/poly4_discr.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Polynomial discrimination</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 8.6.2, Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX by Joelle Skaf - 10/23/05</span>
<span class="comment">% (a figure is generated)</span>
<span class="comment">%</span>
<span class="comment">% The goal is to find the polynomial of degree 4 on R^n that separates</span>
<span class="comment">% two sets of points {x_1,...,x_N} and {y_1,...,y_N}. We are trying to find</span>
<span class="comment">% the coefficients of an order-4-polynomial P(x) that would satisfy:</span>
<span class="comment">%           minimize    t</span>
<span class="comment">%               s.t.    P(x_i) &lt;= t  for i = 1,...,N</span>
<span class="comment">%                       P(y_i) &gt;= t   for i = 1,...,M</span>

<span class="comment">% Data generation</span>
rand(<span class="string">'state'</span>,0);
N = 100;
M = 120;

<span class="comment">% The points X lie within a circle of radius 0.9, with a wedge of points</span>
<span class="comment">% near [1.1,0] removed. The points Y lie outside a circle of radius 1.1,</span>
<span class="comment">% with a wedge of points near [1.1,0] added. The wedges are precisely what</span>
<span class="comment">% makes the separation difficult and interesting.</span>
X = 2 * rand(2,N) - 1;
X = X * diag(0.9*rand(1,N)./sqrt(sum(X.^2)));
Y = 2 * rand(2,M) - 1;
Y = Y * diag((1.1+rand(1,M))./sqrt(sum(Y.^2)));
d = sqrt(sum((X-[1.1;0]*ones(1,N)).^2));
Y = [ Y, X(:,d&lt;0.9) ];
X = X(:,d&gt;1);
N = size(X,2);
M = size(Y,2);

<span class="comment">% Construct Vandermonde-style monomial matrices</span>
p1   = [0,0,1,0,1,2,0,1,2,3,0,1,2,3,4]';
p2   = [0,1,1,2,2,2,3,3,3,3,4,4,4,4,4]'-p1;
np   = length(p1);
op   = ones(np,1);
monX = X(op,:) .^ p1(:,ones(1,N)) .* X(2*op,:) .^ p2(:,ones(1,N));
monY = Y(op,:) .^ p1(:,ones(1,M)) .* Y(2*op,:) .^ p2(:,ones(1,M));

<span class="comment">% Solution via CVX</span>
fprintf(1,<span class="string">'Finding the optimal polynomial of order 4 that separates the 2 classes...'</span>);

cvx_begin
    variables <span class="string">a(np)</span> <span class="string">t(1)</span>
    minimize ( t )
    a'*monX &lt;= t;
    a'*monY &gt;= -t;
    <span class="comment">% For normalization purposes only</span>
    norm(a) &lt;= 1;
cvx_end

fprintf(1,<span class="string">'Done! \n'</span>);

<span class="comment">% Displaying results</span>
nopts = 2000;
angles = linspace(0,2*pi,nopts);
cont = zeros(2,nopts);
<span class="keyword">for</span> i=1:nopts
   v = [cos(angles(i)); sin(angles(i))];
   l = 0;  u = 1;
   <span class="keyword">while</span> ( u - l &gt; 1e-3 )
      s = (u+l)/2;
      x = s * v;
      <span class="keyword">if</span> a' * ( x(op,:) .^ p1 .* x(2*op) .^ p2 ) &gt; 0,
          u = s;
      <span class="keyword">else</span>
          l = s;
      <span class="keyword">end</span>
   <span class="keyword">end</span>;
   s = (u+l)/2;
   cont(:,i) = s*v;
<span class="keyword">end</span>;

graph = plot(X(1,:),X(2,:),<span class="string">'o'</span>, Y(1,:), Y(2,:),<span class="string">'o'</span>, cont(1,:), cont(2,:), <span class="string">'-'</span>);
set(graph(2),<span class="string">'MarkerFaceColor'</span>,[0 0.5 0]);
title(<span class="string">'Optimal order-4 polynomial that separates the 2 classes'</span>)
<span class="comment">% print -deps min-deg-discr.eps</span>

<span class="comment">%%%% Dual infeasible ?????</span>
</pre>
<a id="output"></a>
<pre class="codeoutput">
Finding the optimal polynomial of order 4 that separates the 2 classes... 
Calling sedumi: 228 variables, 17 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 17, order n = 215, dim = 229, blocks = 2
nnz(A) = 3393 + 0, nnz(ADA) = 287, nnz(L) = 152
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            5.60E+02 0.000
  1 :  -4.97E-01 2.36E+02 0.000 0.4220 0.9000 0.9000   2.73  1  1  1.3E+02
  2 :  -1.84E-01 1.51E+02 0.000 0.6389 0.9000 0.9000   5.08  1  1  2.9E+01
  3 :  -2.46E-02 7.59E+01 0.000 0.5024 0.9000 0.9000   4.64  1  1  5.1E+00
  4 :  -1.13E-02 3.19E+01 0.000 0.4205 0.9000 0.9000   1.51  1  1  1.9E+00
  5 :  -6.08E-03 1.48E+01 0.000 0.4643 0.9000 0.9000   0.80  1  1  1.1E+00
  6 :  -8.09E-04 4.48E+00 0.000 0.3024 0.9000 0.9000   0.10  1  1  9.3E-01
  7 :   7.49E-03 2.84E+00 0.000 0.6334 0.9000 0.9000  -0.40  1  1  7.1E-01
  8 :   2.80E-02 1.50E+00 0.000 0.5302 0.9000 0.9000   0.36  1  1  3.8E-01
  9 :   3.40E-02 5.24E-01 0.000 0.3483 0.9000 0.9000   0.52  1  1  1.5E-01
 10 :   3.73E-02 2.14E-01 0.000 0.4080 0.9000 0.9000   0.67  1  1  7.4E-02
 11 :   3.89E-02 1.01E-01 0.000 0.4702 0.9000 0.9000   0.68  1  1  3.9E-02
 12 :   3.98E-02 5.69E-02 0.000 0.5662 0.9000 0.9000   0.60  1  1  2.6E-02
 13 :   4.03E-02 3.12E-02 0.000 0.5481 0.9000 0.9000   0.67  1  1  1.5E-02
 14 :   4.07E-02 1.69E-02 0.000 0.5430 0.9000 0.9000   0.65  1  1  9.2E-03
 15 :   4.09E-02 4.50E-03 0.000 0.2655 0.9000 0.9000   0.86  1  1  2.5E-03
 16 :   4.10E-02 1.14E-03 0.000 0.2530 0.9000 0.9000   0.88  1  1  6.7E-04
 17 :   4.10E-02 2.30E-05 0.000 0.0202 0.9900 0.9900   0.99  1  1  1.4E-05
 18 :   4.10E-02 6.03E-06 0.000 0.2622 0.9000 0.8743   1.00  1  1  3.3E-06
 19 :   4.10E-02 1.39E-06 0.000 0.2306 0.9003 0.9000   1.00  1  1  7.6E-07
 20 :   4.10E-02 1.93E-07 0.000 0.1390 0.9107 0.9000   1.00  1  1  1.3E-07
 21 :   4.10E-02 3.64E-08 0.000 0.1883 0.9106 0.9000   1.00  1  1  3.1E-08
 22 :   4.10E-02 9.83E-09 0.046 0.2700 0.9169 0.9000   1.00  2  2  1.0E-08

iter seconds digits       c*x               b*y
 22      0.1   7.9  4.0993225478e-02  4.0993225002e-02
|Ax-b| =   9.0e-09, [Ay-c]_+ =   0.0E+00, |x|=  5.6e-01, |y|=  1.0e+00

Detailed timing (sec)
   Pre          IPM          Post
2.000E-02    1.300E-01    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 26.9346.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -0.0409932
Done! 
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="poly4_discr__01.png" alt=""> 
</div>
</div>
</body>
</html>